m 

o 

< 



X 



Exact Fourier Spectrum Recovery 

M. Andrecut 
April 9, 2013 

ISIS, University of Calgary, Alberta, T2N 1N4, Canada 
(email: mandrecu@ucalgary.ca) 

Abstract 



r~^ I Discrete Fourier Transform (DFT) is widely used in signal processing to analyze the frequencies in 

a discrete signal. However, DFT fails to recover the exact Fourier spectrum, when the signal contains 
(~| . frequencies that do not correspond to the sampling grid. Here, we present an exact Fourier spectrum 

recovery method and we provide an implementation algorithm. Also, we show numerically that the 
Cw ' proposed method is robust to noise perturbations. 

C^ ' 
^.. 

c/3 ! 1 Introduction 

o 

^^' Discrete Fourier Transform (DFT) is one of the most important tools in digital signal processing. DFT is 

<~i ■ used in both analysis and synthesis of discrete signals and systems, and has frequent applications in optics, 

'^- spectroscopy, magnetic resonance, quantum computing, seismology and astrophysics [IJE]- In particular, 

DFT allows discrete-time signals and systems to be analyzed in the frequency domain. 
^ I DFT is an invertible linear operator T : C^ — > C^, that transforms a set of complex numbers, defining 

^V I the signal z„ — z(i„) G C , sampled on the time grid t„ = nAt, n = 0, ..., N — 1, into another complex set 

f^ ' of numbers: 

(N' 

^' 

o 

CO ' where Vm — m/{NAt), m — 0, ..., N — 1 is the frequency sampling grid. The inverse DFT is defined as 

following: 



N-l 
ra=0 



N-l 



Zn — ^n [{Zm}m=a\ ~ /_^ ^m^ " " . (2) 

m— 

^.' A major short coming of DFT is that it correctly recovers the Fourier spectrum only if the signal contains 

frequencies from the sampling grid {Dm — i^/(.^^t)}m=o HI IS- ^^ Figure 1 we give a DFT failing recovery 
example where the signal has a Fourier spectrum containing K — 8 components with frequencies Vk {k = 
0, ...,K — 1), that are not included in the sampling grid, i.e. v^ ^ {vm = m/{NAt)}^-^'lQ. One can see that 
the spectrum obtained using DFT does not even resemble the original frequency spectrum embedded in the 
signal. The standard approach to overcome this problem is to sample the signal at a larger number of points 
N , such that the frequencies in the sampling grid begin to approximate well the frequencies embedded in 
the signal. This approach is obviously not feasible if one can sample the signal only at a fixed number of 
points, constrained by the experimental setting for example. The main question we would like to address 
here is how to recover exactly the Fourier spectrum if the signal in question contains frequencies that do not 
fall on the sampling grid? Here, we derive an iterative method which recovers exactly the Fourier spectrum 
embedded in the signal, and we provide an implementation algorithm which is robust to noise perturbations. 



2 Exact Recovery Method 

The goal of the Fourier spectrum recovery problem is to exactly find the Fourier components {Z^, '^fe};;^ , 
K < N/2, from the signal: 

K-l 

zn = Y. Zke^^''""" , n = 0, ..., iV - 1, (3) 

k=0 

including the case when the frequencies Vk do not fall on the sampling grid. 

We consider an iterative method, and we assume that initially Zk — 0, Vk — 0, k — 0,...,K — 1, 
and the residual is f = [zq, ..., zn-i] G C^. At every step the method uses a nonlinear optimization 
strategy to determine the amplitude Zk and the frequency Vk corresponding to the Fourier component 
^{vk) = [e^'^**"'^*, ..., e^'^**™"^'^''] e C^, which minimize the residual in the least squares sense. A cyclic 
strategy is also employed to reassess and correct the {Zk, Vk) parameters, while holding the others fixed. In 
order to reassess the component k with the parameters {Zk, Vk) we remove it from the solution, by adding 
its contribution to the residual: 

f ^ f + Zki){vk)- (4) 

Here, we have introduced the assign operator variable ^ expression / program (re-sets the variable with 
the returned value of the expression / program) . Now, we compute another set of parameters {Zk, Vk) which 
minimize the residual f: 

(Zfe,i/fc) ^ argminS'fc(Z, i^), (5) 

where 

2 ^~1 

Sk{Z,v) = ||r - Zi:{y)\ - ^ [r„ - Ze^-*""-]' . (6) 

n=0 

By solving the above minimization problem we find a new set of parameters Zk and Vk , which will produce 
the residual: 

f ^f - Zk-ip{vk)- (7) 

The new residual f will be smaller or at least equal with the previous one, due to the intrinsic convergence 
mechanism of the method (as shown below). The worst case would be when the values of the reassessed 
parameters {Zk,Vk) are identical with the initial ones and no correction can be made. Obviously, this 
procedure can be repeated cyclically for each component until the norm of the residual becomes zero, or 
smaller than some prescribed positive threshold: ||f|l < £,, ^ 0. When the signal is contaminated with 
Gaussian noise with zero mean and standard deviation cr, the stopping threshold should be comparable with 
the standard deviation of the noise, i.e. e^ ~ a. 

The major difficulty of the method is solving the successive nonlinear minimization problems (5). The 
functions Sk{Z,v) are highly nonlinear, with complex arguments and multiple local minima, which makes 
them very difficult for the existing algorithms [5]. Here, we develop a new iterative algorithm to solve the 
nonlinear minimization problem. First, we need to find an initial estimate Vk- This is is done simply by 
solving the problem: 



Vk ■^ argniax 



f, 1p{Vm) 



(8) 



where (., .) is the standard inner product operator in the complex Hilbert space, and Dm are the frequencies 
from the sampling grid {P,„ — m/{NAt)}^'Zl. We also notice that the amplitude Zk and the frequency 
i>k can be "decoupled", in the sense that for a given Vk, one can easily calculate an estimate of Zk as the 
projection of f on ^{vk)'. 

Zk=N-Uf,^{vk)). (9) 



This estimate of Z^ guarantees the decrease of the residual (7), since we have: 



r-Zktpii^k) = (f - Zki){vk),f ~ Zkil){vk) 



{r,r) - Zl (^f,ip{iyk)J - Zk{i}{vk),fj + \Zk\ {i){vk)A{vk) 
\\r\\^-N\Zk\^-N\Zk\'+N\Zk\^<\\f\\\ 



Reciprocahy, ^{vk) is orthogonal to the residual (7): 

■ip{i^k),f - Zki^ivk) 



(10) 



i'{i'k),rj - Zl (^^p{i^k),ip{vk] 
N4{t)^NcUt)^0. 



(11) 



Now, let us assume that Avk is the unknown correction to Vk in the next iteration step. Therefore, we 
have: 



with 



and from here we easily find: 



^p{i^k + Az^fc) ~ ^p{i^k) H 1 Ai^fe, 

dv 



dVyVk) n^-\, 2Trw^ta I „27rij/ttAr-ll^ 

— == iTTi [toe *= ", ...,i7v-ie <- « ij , 



Avk 



Zk 



dip{vk) 



dv 



f - Zk'4){vk),Zk 



dip{vk) 
dv 



Again, this estimate of the correction l\vk guarantees a further decrease of the residual: 

r ^r - Zkip{Vk) - Zk — CxVk, 

dv 



since we have: 



r - Zkip(Vkj - Zk — Az/fe 

dv 



r - Zk-4>{vk) 



-Avl 



dipivk) 



dv 



< 



f - Zkip{vk) 



Reciprocally, the direction Zk ^y''' is also orthogonal to the next residual: 



di[j(vk) r7 7 / dtp(vk) 

Zk — 1 ,r- Zktp{vk) - Zk — -: Ai/fc 

dv dv 



(12) 
(13) 

(14) 



(15) 



(16) 



Zk ; ,f ~ Zki>{vk) ) ~ Avk ( Zk "'^]' '^' ,Zk'—, / 

dv \ dv dv 



d-ipivk) y dip{vk) \ 



Avk 



d-4;{vk) 



dv 



- Avk 



Z, 



dip{vk) 



dv 



0. 



After /\vk is determined, we update the frequency of the current component using Vk 



(17) 
Vk + Re{/\vk} 



Algorithm 1 Exact Fourier Spectrum Recovery (EFSR) metlrod 



z; signal 

K; number of components to be recovered 

Er] admissible residual level 

e^; admissible frequency error 

S\ maximum number of correction cycles 

T; maximum number of iterations per cycle 

f ■<— z; initial residual 

s ^— 0; cycle index 

do{ 

k ^ mod(s, K); 



Vk ^ arg maxp^ 


('■■ 


^{ym)) 


i^O; 


do{ 


Zk^N-^{r,i^{yk))-. 


Ai^k^ 






-'\^- 



f - Zkij{i^k),Zk-^ 



dipjuk] 



Vk-^ Vk + Re{/\vk}] 

t^t + 1: 

}while(|Ai/fe| > £0 and t < T) 

f ^r- Zktpiiyk); 
s ^ s + 1; 

}while(||f|| > Sr and [s/K\ < S) 
return {Zk^i^kJkJo^ 



(since the frequency must be a real number) and the residual: f ^^ r — Zkip{vk), and we perform a new 
iteration. The iterations stop when [Avkl < Ci/, where < e,y -C 1 is a small acceptable error, or when the 
number of iterations exceed a maximum number. Therefore, at every step, the method finds two directions 
ipiiyk) and respectively Zk-^^^j^, which are orthogonal to the next residual. The parameters (Zkji'k) are 
the projection of the current residual on these two directions. Thus, as a consequence of the above proof, 
the method is guaranteed to converge. The pseudo-code of the Exact Fourier Spectrum Recovery (EFSR) 
algorithm is given in Algorithm 1. 

In the example form Figure 1 one can see that DFT cannot recover the spectrum at all, while the EFSR 
method recovers exactly all the components in the Fourier spectrum (e,y = Sr = 10^®, T = 10^, S = K"^) 
in about s = 10 iteration cycles. We should mention that the EFSR method will recover exactly all the 
cases when the frequencies from the spectrum are well separated among them by at least the Nyquist limit, 
which is the highest frequency that can be coded at a given sampling rate At in order to be able to fully 
reconstruct the signal, i.e. Av ~ 2.0/ (N At). Thus, for a successful recovery we must have: {I'k ~ ^j\ ^ ^^) 
for k^j, k,j ^0,1,...,K -1, K < N/2. 

In order to investigate the effect of the noise perturbation we consider the perturbed signal Zn + rjn, where 
r]n is the measurement noise, i.e. a random variable normal distributed with mean (?7„) = and standard 
deviation a. Also, as a measure of noise contamination we use the signal to noise ratio defined as: 



SNR- 



fRMS. 



signal 



\RMS„ 



a RMSgign^i, 



(18) 



where RMS is the root mean square. 

In Figure 2 we give an example with the signal to noise ratio: SNR = 10. One can see that the EFSR 
method still recovers the spectrum almost exactly. In this case the iterations were stopped when the norm of 
the residual becomes comparable with the standard deviation a of the noise in the measured data, i.e. when: 



||f|| < a, and thus e^ = cr. In order to estimate the recovery capabihties of the EFSR method in the presence 
of noise, we define the foUowing average relative errors between the recovered spectrum {Z^ , v)[ }y.^i^ and 



the original spectrum {Zk,i'k}k^Q'- 



£^ = 17 E ^^ X 100%. (19) 
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K-l 

^- 

k=0 


4^^- 


Zk 


K 


\Zk 




1 


K-l 

>: 

fc=0 


4' - 


Vk 


K 


i^k 





X 100%. (20) 

The recovery errors as a function of the SNR are represented in Figure 3. The computation was performed 
by averaging over M = 10^ cases for each SNR value, using the following parameters: N = 128, -ftT = 8, 
T = 10^, S = K^ . The numerical results show that the method is robust to noise perturbations, and recovers 
the spectrum well up to SNR ~ 10, when begins to deteriorate. Also, it is interesting to note that the error 
of the complex amplitudes ez grow much faster than the frequency error Su. This is due to the fact that by 
increasing the noise perturbation, the information about the phase of the amplitude is gradually lost, while 
the frequencies are still well recovered even for up to SNR = 2. 

3 Application to Faraday Rotation Measure Synthesis 

Faraday rotation is a physical phenomenon where the position angle of linearly polarized radiation propagat- 
ing through a magneto-ionic medium is rotated as a function of frequency. Recently, the Faraday rotation 
measure (RM) synthesis has been re-introduced as an important method for analyzing multichannel polar- 
ized astrophysical radio data, where multiple emitting regions are present along the single line of sight of 
the observations. In practice, the method requires the recovery of the Faraday depth function from mea- 
surements restricted to limited wavelength ranges, which is an ill-conditioned deconvolution problem, raising 
important computational difficulties (see |U[S] for a detailed description). 

Faraday rotation is characterized by the Faraday depth (in radm~^), which is defined as: 

f^observer 

(j){r) = 0.81 / n^B ■ dr, (21) 



source 



where rie is the electron density (in cm '^) , _B is the magnetic field (in ^G), and dr is the infinitesimal path 
length (in parsecs). We also define the complex polarization as: 

P{X')^Q{X'')+iU{X^), (22) 

where /, Q, U are the observed Stokes parameters. The observed polarization P(A^) originates from the 
emission at all possible values of the Faraday depth 0, such that: 

P(A2) = / F(0)e2'^^ #, (23) 

where F{(j)) is the complex Faraday depth function (the intrinsic polarized flux, as a function of the Faraday 
depth). Thus, in principle F((/)) is the inverse Fourier transform of the observed quantity P(A^): 

P((/.) = / P{\^)e-'^"f'^ d\\ (24) 



In general, the number of observations N is limited by the number of independent measurement channels, 
and therefore there are many different potential Faraday depth functions consistent with the measurements. 
The usual approach to resolving such ambiguities, is to impose some extra constraints on the Faraday depth 
function. Here, we consider that the complex Faraday depth function is approximated by a small number of 
components. More specifically, we assume that F((j)) contains K (unknown) Dirac components fk5{4> ~ 4'k)i 
with complex amplitudes f^ G C, and depths 0^ G M, A; = 0, ...,K — 1: 

K~l 



k=0 



K-l K-l 



F{(f)=Y.h5{^-<l^k), (25) 

Taking the Fourier transform of F{(l)) we obtain: 

E fk^i^ - Me'^^^^dc^ = Y. he'"^"^^- (26) 

-°° fc=0 fc=0 

The goal of the RM synthesis is to find F(0) from the observed values P{X^) = Pn (i.e. Qn and [/„) over N 
discrete channels A^ G ['^mim^max]^ n = 0,1,...,N — 1. From the numerical point of view one can apply the 
EFSR method directly, by substituting the Faraday depth and the squared wavelength with: O tt;^ and 

In Figure 4 we give such an example simulated for the Westerbork Synthesis Radio Telescope (WSRT) 
experiment layout P| fS] . The various parameters associated with the WSRT experiment layout are listed 
bellow: 



Wave length range: A^„ = 0.639 m^, A^^^^ = 0.905 m^; 

The width of an observing channel: 6X^ ~ 0.002f m^; 

Maximum observable Faraday depth: cpmax — SOOradm^^; 

Depth space resolution (equivalent to half Nyquist frequency): Scj) ~ fSradm"^; 

We randomly generated K ^ 5 sources with complex amplitudes, Faraday depth —(j)max < 4> < (pmax- 
All the components are separated at the Nyquist limit. We also added noise to the generated polarization 
vector P(A^), such that SNR = fO. One can see that while DFT fails to recover the depth spectrum, while 
the EFSR method correctly determines the all the components embedded in the polarization signal, with a 
small error in the phase, due to the ambiguity induced by the noise in the data. 

4 Conclusion 

We have presented an exact Fourier spectrum recovery method, which can be applied in all cases, including 
those in which the signal contains frequencies that are not corresponding to the sampling grid. Also, we have 
provided an implementation algorithm and we have shown numerically that the method is also extremely 
robust to noise perturbations. 
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Figure 1: An example of noiseless signal {SNR — oo), sampled at A^ = 64 time steps (1st line), containing 
K — 8 Fourier components with frequencies not falling on the sampling grid (2nd line). The DFT method 
is failing in this case (3rd line), while the EFSR method recovers the Fourier spectrum exactly (4th line). 
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Figure 2: An example of noisy signal {SNR — 10), sampled at A^ = 64 time steps (1st line), containing 
K = 8 Fourier components with frequencies not falling on the sampling grid (2nd line). The DFT method 
is failing in this case (3rd line), while the EFSR method recovers the Fourier spectrum almost exactly (4th 
line) . 




Figure 3: Recovery errors as a function of the SNR: the error of the complex amphtudes ez grow much 
faster than the frequency error e^. 
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Figure 4: An example of noisy rotation measure synthesis {SNR — 10), for the WSRT experiment layout 
(1st line), containing K — 5 components (2nd line). The DFT method is failing in this case (3rd line), while 
the EFSR method recovers the spectrum almost exactly (4th line). 
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